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Abstract 


This paper investigates the effectiveness of employing the traditional Leverett approach for describing the capillary-induced flow in thin-film fuel 
cell diffusion media (DM) with mixed wettability. A one-dimensional steady-state analytical model is developed to analyze the capillary transport 
of liquid water through the thin-film DM. A capillary pressure—liquid saturation relation is derived based on the experimental data reported by 
Gostick et al. [J.T. Gostick, M.W. Fowler, M.A. Ioannidis, M.D. Pritzker, Y.M. Volfkovich, A. Sakars, J. Power Sources, 156 (2006) 375]. The 
empirical fit is applicable to the fuel cell DMs tested by reference [J.T. Gostick, M.W. Fowler, M.A. Ioannidis, M.D. Pritzker, Y.M. Volfkovich, 
A. Sakars, J. Power Sources, 156 (2006) 375] under a limited set of conditions, providing a means to evaluate the effectiveness of the traditional 
Leverett approach. Furthermore, a new characteristic relative permeability correlation appropriate for the tested DMs under the given conditions 
was obtained by fitting the experimental capillary pressure data [J.T. Gostick, M.W. Fowler, M.A. Ioannidis, M.D. Pritzker, Y.M. Volfkovich, A. 
Sakars, J. Power Sources, 156 (2006) 375] into four well-established empirical models. The empirically-derived relationships are then integrated 
into an analytical model framework in order to compare the liquid saturation profiles predicted by both approaches. The results show that use of the 
standard Leverett approach equipped with Leverett J-function can lead to significant errors; therefore, an extension of this approach appropriate 
for fuel cell DM with mixed wettability is needed for reliable model predictions. Finally, a sensitivity analysis is performed to assess the relative 


significance of the various input parameters on the predicted saturation profiles. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The efficient operation of polymer electrolyte fuel cells 
(PEFCs) requires complex water management strategies [1-5] 
due to the interaction of electrochemistry, material science, and 
transport phenomena. Reliable and efficient operation requires 
a delicate balance between membrane hydration and the avoid- 
ance of cathode flooding. Adequate water in the membrane 
electrolyte and catalyst layers is important for achieving suit- 
able performance. However, the presence of excess liquid water 
in the diffusion media (DM), catalyst layer and the flow channel 
increases the reactant transport resistance to the electrode cata- 


Abbreviations: DM, diffusion media; J(s), Leverett function; PEFC, polymer 
electrolyte fuel cell; MPL, micro-porous layer 
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lyst surface, causing severe performance loss. The existence of 
two-phase flow in the porous DM presents a major bottleneck in 
the design optimization of PEFCs due to the highly interrelated 
and ill-defined transport phenomena involved. 

Because of the minute length scales and presumed anisotropic 
nature of the thin-film DM, the multi-phase transport behav- 
ior inside the DM is not yet fully understood. Within the DM, 
gravitational and viscous forces are relatively small compared 
to capillary forces [5]. Hence, two-phase transport in the DM is 
primarily controlled by capillary action. The pressure difference 
arising from the interfacial tensions at the liquid—gas interface 
front is referred to as the capillary pressure. The capillary pres- 
sure is a strong function of liquid saturation (sı) which is a key 
parameter representing the available pore space, through which 
reactant gases can diffuse. 

Fuel cell performance is known to exhibit hysteresis, for a 
variety of reasons, but often because of liquid water motion. 
Many studies [5—9] reveal that liquid water transport is a dom- 
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Nomenclature 


Faraday constant (96,485 C (eq mol)~!) 
current density (A cm~?) 
permeability (m7?) 
relative permeability (unitless) 
fitting parameter (unitless) 

W molecular weight (g mol~!) 
pressure (Pa) 
liquid saturation (unitless) 
phase velocity in x direction (m s7!) 
distance, m 
phase velocity (m s7!) 
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Greek letters 


a net water transfer coefficient (unitless) 
E porosity (unitless) 

y surface tension of water—air (N m7!) 
Xr pore distribution index (unitless) 

u absolute viscosity (kg m s7!) 

v kinematic viscosity of water (m? s7!) 
0 contact angle (°) 

p density (kg m7?) 

Subscripts 

C capillary 

e effective 

8 gas 

int interface 

irr irreducible 

l liquid 

w wetting phase 

nw non-wetting phase 


inating factor influencing the performance of the fuel cell, and 
is closely related to the liquid saturation distribution in DM. 
Therefore, accurate prediction of the liquid water distribution in 
the DM as a function of material properties is of great interest 
to understand the capillary transport phenomenon in the DM 
[6-11]. An extensive review of the recent modeling efforts is 
provided in Ref. [12]. Since the vast majority of the background 
work on defining the capillary pressure—saturation relation has 
been performed in the soil science and oil extraction fields, a 
generic Leverett approach from soil science has been commonly 
employed to represent the water transport or retention behavior 
of the fuel cell DM as a first step toward achieving an accurate 
two-phase transport model for fuel cells studies [5-7]. The Lev- 
erett approach, indeed, is viewed as an indispensable tool for 
many modelers and serves as a useful starting point to model 
the liquid transport in the porous diffusion media. Although the 
usefulness of the Leverett approach is clear, the degree of the 
applicability of the traditional Leverett approach in its original 
form to the highly anisotropic thin-film DM has not been con- 
clusively established. Since the naturally hydrophilic fuel cell 
DM material itself is impregnated with an anisotropic coating 


of hydrophobic material (PTFE or other), it yields a complex 
bi-modal pore size distribution with mixed wettability, compli- 
cating the transport phenomena therein. 

Specific concerns center around the characteristic morpho- 
logical differences in fuel cell DM and conditions under which 
the Leverett approach equipped with J-function was derived, 
including the fact that it was: empirically derived to simulate 
a range of common soils with uniform wettability, ignoring 
the observed hysteresis in porous media and most importantly, 
developed for isotropic soil beds with a high volume to surface 
area ratio [13]. As a result of the anisotropic coating of hydropho- 
bic material (PTFE or other), the liquid water can be both the 
wetting (on untreated DM fiber), and non-wetting phase (on 
PTFE coated fiber) phase [12], complicating the mathematical 
description of the transport phenomena. Furthermore, the defini- 
tion of contact angle representing the wettability characteristics 
of the anisotropic DM with mixed wettability is typically taken 
as a Statistical average over the entire medium, obscuring local 
effects which may differ from the whole. Therefore, the appli- 
cability of the traditional Leverett approach to thin-film fuel cell 
DM is questionable, and is the focus of this paper. 

Another important parameter governing the capillary trans- 
port inside the DM is relative permeability, which is also a 
strong function of liquid saturation. In PEFC modeling studies, 
a typical relative permeability expression for non-consolidated 
sands is adapted to describe the phase transport in DM due to its 
simplicity, that is: 


Ley (1) 


where kr-nw and Spw represent the relative permeability of 
non-wetting phase and non-wetting saturation, respectively. 
However, the applicability of this correlation to the thin-film 
DM is also questionable and yet unproven. Direct experimen- 
tal evidence is still needed to verify the applicability of these 
empirical correlations, since the reliability of the computational 
output is directly related to the quality of the input parameters. 

The direct measurement of liquid saturation distribution 
through the pores of the DM is an experimental challenge and 
only limited experimental studies are yet available in the open 
literature [1,2]. Recently, Gostick et al. [10] reported a novel 
experimental approach for determining the dependence of cap- 
illary pressure on wetting phase saturation along a desaturation 
path (drainage) for various DMs. This technique, the method 
of standard porosimetry (MSP), is capable of distinguishing 
the hydrophilic pore distribution from the total pore network, 
unlike mercury intrusion porosimetry. Using standard porosime- 
try, Gostick et al. [10] measured the static capillary pressure as 
a function of saturation during the drainage process for a certain 
PEFC DM materials under a limited set of conditions. 

The present study reports a detailed comparative analysis 
of the constitutive relations! relations describing the capillary- 
induced flow in fuel cell DM, particularly focusing on the 
applicability of the commonly employed standard Leverett 


' The term constitutive relation is used here to denote the empirical relation 
of capillary pressure and liquid saturation in a fuel cell diffusion media. 
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Fig. 1. Schematic of control volume. 


approach to thin-film fuel cell DM with mixed wettability. The 
saturation prediction of the Leverett approach in its original form 
(as applied to the most fuel cell models) was compared with the 
prediction of the empirical capillary pressure—saturation corre- 
lation derived from the experimental data provided by Ref. [10]. 
The reasons underlying the variations in liquid saturation pre- 
dictions within the DM were linked to the prediction potential 
of the Leverett approach. With the aid of these analyses, the 
ineffectiveness of employing Leverett approach in its original 
form is discussed and the necessity of development of a modi- 
fied Leverett approach appropriate for thin-film fuel cell media 
is addressed. 


2. Theoretical 
2.1. Theoretical analysis of liquid transport through DM 


A basic analytical model describing the capillary liquid flow 
has been developed to predict the liquid saturation distribution 
inside the DM. This model is constructed to provide a tool to 
assess the saturation predictions of the constitutive relations 
while isolating all other effects, including: humidity, 2D effects, 
inlet cell conditions and time dependent parameters, etc. Fig. 1 
shows a schematic of the chosen control volume in the fuel cell. 
In the present model, the liquid water flowing inside the DM is 
assumed to be Newtonian and incompressible, and the flow in 
the pores of the DM is presumed to be uni-directional, steady 
and laminar. 

In the steady state, once the gas phase is fully saturated with 
water vapor, assuming evaporation and condensation do not take 
place in the DM, liquid water flow becomes the only mode of 


1.417(1 — sy) — 2.12001 — s1) + 1.263(1 — s1)? 
J(sı) = 


1.4175) — 2.1205? + 1.2635; 


the water transport across the DM [6]. Water generation in the 
catalyst layer, water transport across the membrane, and conden- 
sation or evaporation within the DM pores cause non-uniformity 
in the saturation distribution. As hydrophobic DM pores are 
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Fig. 2. Illustration of capillary-induced liquid flow in a fuel cell DM. 


filled by the liquid water, the liquid phase pressure increases, 
eventually driving the liquid water from higher to lower liquid 
pressure regions. As aresult, the liquid water in the DM is driven 
via capillary action produced by the liquid saturation gradients 
in DM. This phenomenon is illustrated in Fig. 2. 

In a hydrophobic DM, capillary pressure is defined as the 
difference between liquid pressure (non-wetting phase) and 
reactant gas pressure (wetting phase). 


Po = P — Pg (2) 


The liquid flow through the porous DM can be described by 
Darcy’s law. Darcy’s law relates the flow rate through a porous 
medium to the pressure gradient across the porous media and the 
permeability of the porous media. Assuming constant gas phase 
pressure (Pg © const) and no convective gas flow, the governing 
liquid transport equation in the porous DM can be expressed as: 


3 —k k; > > > 
v = —VPc_ where VPe ~ VP, (3) 
u 


where k is absolute permeability of the porous media, kr the fluid 
relative permeability, u the absolute viscosity of the fluid and 0, 
is the velocity of the flowing fluid. 

A generic Leverett approach has been commonly employed to 
describe the capillary pressure—saturation relation of the DM in 
two-phase PEFC models. Leverett [21] related capillary pressure 
to saturation in porous media by: 


1/2 
Po=y cosø( £) Is) (4) 


where k, £, 0, y and J(s1) represent the gas phase or liquid abso- 
lute permeability, porosity of the DM, contact angle between 
liquid—solid phase, surface tension and the Leverett function in 
terms of liquid saturation, respectively. The standard Leverett 
function, J(s}), is given as [21]: 

if 0 < 90° => Hydrophilic 6) 
if 6 > 90° = Hydrophobic 


The liquid pressure (P}) is a function of local liquid saturation; 
hence, the capillary pressure gradient can be expressed as: 


> > e /(dJ\ = 
VP ~ VPc=y cosh,/— | — S| (6) 
k ds, 
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Inserting Eq. (6) into Eq. (3), the governing equation yields: 


F kk, 9 e /dJ\ = (7) 
= cos, / — | — 
“l W Y k \ds È 


Assuming all the water produced by electrochemical reaction 
in a PEFC is in liquid form, the mass flux of liquid water going 
into the cathode DM is equal to sum of the net water transport 
from the anode side and the amount of water generated in the 
catalyst layer due to the electro-chemical reaction. The net water 
flux across the membrane from anode to cathode can be charac- 
terized by a net water transport coefficient, œ [29]. Combining 
the net water transport with the water production, the net water 
flux to the cathode DM yields: 


I 
pv, = apm + MWO (8) 


where MW!2° is the molecular weight of the water, œ the 
net water transport coefficient to the cathode DM, Z the cur- 
rent density and F is the Faradic constant (96,485 C mol eq~'). 
Assuming one-dimensional flow, substituting Eq. (8) into Eq. 
(7), the final form of the governing equation becomes: 


T2a+ VMWPOn, dJ ds 0) 
2F yp cos0v ke 7 "dsi dx 


This relationship describes the steady-state liquid water 
transport across the fuel cell DM. As seen from Eq. (9), the 
operational conditions and DM material properties are coupled 
with the fluid properties and the mass transport rate are explicitly 
linked to the liquid saturation. 


2.1.1. Micro-porous layer interface 

Generally, the cathode side of PEFC consists of a coarse 
DM and a finer micro-porous layer (MPL), both of which are 
generally classified as hydrophobic, based on surface contact 
angle. However, the DM is actually a mixed wettability network, 
since the base material of carbon-fiber is hydrophilic. Micro- 
porous layers are introduced between the DM and catalyst layer 
(CL) to provide better membrane electrical contact and improve 
the water management. The liquid water transport of the MPL 
is governed by capillary action. In equilibrium, the capillary 
pressure across the interface between the DM and MPL is con- 
tinuous. The different material properties of these two layers 
cause a discontinuity in the liquid saturation across the interface 
[7], as illustrated in Fig. 3. The magnitude of the discontinuity 
or the jump in the saturation at the interface strongly depends 
on the material properties of these two layers, specifically the 
hydrophobicity and pore radius. Imposing the capillary pressure 
balance at the interface yields: 


DM|__ pMPL 
P C |= Pe ed Interface 


1/2 1/2 

EMPL EDM 

COSÊMPL (z=) JME) = cos6pm ( ) IA 
kMPL kpM 


(10) 


With a known liquid saturation at the DM interface, the equi- 
librium MPL liquid saturation at the interface can be calculated 


=P 
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at the interface 
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Fig. 3. Boundary conditions and saturation jump at the interface of the DM and 
MPL. 


from Eq. (10). Since the structural characteristics of MPL and 
physics of the flow through the MPL are similar to the DM, the 
governing equation (Eq. (9)) can also be applied to the MPL in 
order to map the saturation distribution profile inside the MPL. 


2.2. Constitutive equations 


2.2.1. Relative permeability 

When two or more fluids are present in the porous media, 
an additional parameter, known as relative permeability, is used 
to represent the ability of the porous media to conduct various 
fluids simultaneously. Relative permeability is defined as the 
ratio of the effective permeability of one fluid to the absolute 
permeability of the porous media [13,14]. Direct experimen- 
tal measurement of relative permeability for different types of 
soils has been performed. Mathematical approximations based 
on previous experiments are generally preferred for estimating 
the relative permeability because of the difficulty in conducting 
direct experiments. 

Various researchers have proposed correlations based on 
mathematically-derived relationships or experimental data to 
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predict the relative permeability [15]. Most of the existing rela- 
tive permeability correlations are based on the following types 
of physical models: capillary model, statistical model, empiri- 
cal model and network model [15]. Empirically-derived models 
based on direct experimental observation provide the most rea- 
sonable approximations due to their simplicity and reliance on 
direct measurements compared to other types of models [15]. 
Even though all of the above-mentioned models have been used 
to predict the relative permeability [16], they may not be appro- 
priate for use in a fuel cell DM. Typically, the general shape of 
the relative permeability curves is estimated by the following 
equations [15]: 


Kyaw = A(Snw)" 


where A, B, n and m are constants depending upon the structure 
of the porous media. 

In the present study, four well-established empirical correla- 
tions describing the water—air two-phase system (Wyllie model, 
Corey model, Brooks-Corey model and the Van-Genuchten 
model) are used to predict a representative non-wetting phase 
(liquid) relative permeability for differently engineered DMs, 
based on the experimental capillary pressure data provided by 
Ref. [10]. The relative permeability calculation in each model 
is based on the effective saturation (se), which requires the 
irreducible water saturation (Sirr) as an input. The irreducible 
saturation or immobile saturation, (Sirr), is the amount of trapped 
water in the pores of the DM. It can be considered as a threshold 
point below which liquid remains immobile. 


and kw = B(1 — Saw)” (11) 


= = (12) 
1- Sirr 


2.2.2. Wyllie model 

Wyllie [14,15] suggested simple liquid—gas relative perme- 
ability equations for well and non-consolidated sands, based on 
the capillarity concept developed by Purcell [17] and Burdine 
[18]. The Wyllie approach is widely applied in fuel cell modeling 
studies due to its simplistic application [12]. 


kr-w = (1 — Sey (13) 


where Se, nw and w represents the effective saturation, non- 
wetting phase and wetting phase, respectively. 


kr—nw = Ger 


2.2.3. Corey model 

Purcell developed an analytical expression to compute the 
relative permeability from experimental capillary pressure data 
[16]. Burdine [18] introduced the tortuosity factor concept into 
the original capillary pressure model developed by Purcell [17] 


Table 1 
Summary of relative permeability models utilized in this study 


and extended this equation to calculate the multiphase relative 
permeability. Corey combined these two approaches into a sin- 
gle form, based on the assumption that capillary curves over a 
certain range of saturations can be approximated using the linear 
expression, | / P2 = C - Se—w, Where Cis a constant, and obtained 
the following relative permeability equations for non-wetting 
and wetting phases [15]: 


knw = [1 = ew) J = sew]? kw = (sew). (149 
2.2.4. Brooks-Corey model 

Brooks and Corey modified the original Corey capillary 
pressure-saturation relationship to predict the relative perme- 
ability for any pore size distribution [15]. 


knw = (1 — sew) [1 — (sew) >>] 
krw = (sy) AA (15) 


where À is the pore size distribution index of the porous media, 
which is a specific characteristic of the pore size distribution, 
and it is typically determined by fitting the water retention data 
to the Brooks-Corey correlation. 


2.2.5. Van Genuchten model 

One of the most recent relative permeability correlations was 
proposed by Van Genuchten [19,20]. Van Genuchten derived 
an empirical correlation that accounts for variable pore size and 
connectivity concepts of the porous media. This model unites 
the entire range of saturation in a single curve, even though 
not all of the saturation regions (e.g., dry end, wet end) are 
governed by the same physics [19,20]. According to the Van 
Genuchten model, the relative permeability of the non-wetting 
and wetting phases can be calculated as follows: 


Kyaw = (1 — Sew) PL = ewer 

lew = (sE = = Gow" (16) 
where m is a fitting parameter and can be extrapolated from 
capillary pressure—saturation data. The models utilized in the 
present study are summarized in Table |. In all these models, 
the capillary pressure and the relative permeability are coupled; 
thus, each model requires the knowledge of capillary pressure as 
a function of water saturation. However, the standard porosime- 
try technique is incapable of distinguishing the isolated pores 
[10], thus the capillary pressure data presented in Ref. [10] is 
associated with the total saturation. Therefore, a lower limit 
value of Sir, 0.05, is used in the present study, assuming all 
the pores are almost connected, to minimize the isolated pore 


Models 


Relative permeability (non-wetting phase) 


Relative permeability (wetting phase) 


kr—nw = len) 


Wyllie model 
Corey model 


kr-nw = [1 = (Sew) IE = Sew? 


Brooks-Corey model 


Van Genuchten model 


knw = (1 — sew) [L — (Sew)? >] 
kr-nw = (1 — Sew) [1 = (Se-w)!/™"] 


kw = ad = Sei) 
kr—w = (Se—w)4 PRN 
kr—w = (Sexy) 


o kw = (Sew) = (1 = (sew) "T 


Se represents the effective saturation. 
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Table 2 
Material properties of the tested DMs adapted from Ref. [10] 


Material Type Thickness (um) PTFE (wt.%) Hydrophilic porosity Total porosity 
SGL 10BA Paper 380 5 0.63 0.91 
SGL 10BB Paper w/MPL 420 5 0.45 0.87 
Toray 090 Paper 190 0 0.63 0.79 
E-Tek Cloth Cloth 350 0 0.74 0.72 


size effect and estimate the bulk multi-phase connectivity of a 
typical fuel cell DM. 


2.2.6. Empirical fit 

The capillary pressure—saturation curves obtained by Gostick 
et al. [10] were used as benchmark data to obtain a polynomial 
fit explicitly relating the Pc versus sı for the tested fuel cell 
DMs under a limited set of conditions. Four different types of 
diffusion media including: SGL 10BA, SGL 10BB, Toray 090 
and E-TEK Cloth A, were selected, and the material properties 
of these DMs including the total and hydrophilic porosities are 
tabulated in Table 2 (adapted from Ref. [10]). 

Fig. 4a depicts the measurements of capillary 
pressure-saturation (provided by Ref. [10]) for the DM 
chosen in this study. As seen in Fig. 4a, the measured capillary 
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Fig. 4. (a) Measured capillary pressure—saturation data for tested DM sam- 
ples (adapted from Gostick et al. [10]) and (b) comparison of the curve-fit and 
experimental data provided by Ref. [10]. 


pressure for all DM samples follows the same qualitative 
trend, even approaching similar capillary pressures especially 
at low saturations (Spy <0.6). The difference in capillary 
pressure can be attributed to the differences in the specific 
morphological and characteristics and PTFE loading of the 
chosen DMs. As also stated in Ref. [10], a sudden increase 
in capillary pressure is observed after the threshold saturation 
value of 0.8. The sudden change in the measured data at high 
saturation values (Spw>0.8) causes a significant amount of 
scattering, which in turn, complicates the representation by an 
accurate curve-fit. Typically, a fuel cell operates in a saturation 
range of 0-0.6, therefore the capillary pressure measurements 
within the saturation range of 0-0.8 were considered to be 
appropriate benchmark data. A single equation representing 
the actual Pc — sı relation of the tested DMs (Eq. (17)) was 
deduced based on the best fit of the data (0<Spy<0.8) and 
shown in Fig. 4b. Even though the empirical curve-fit starts 
to deviate from the actual data at high saturations (spw >0.8) 
due to the sudden increase in the measured capillary pressure, 
the curve-fit reasonably predicts the measured data within 
the saturation range of 0-0.8 covering the range, in which a 
polymer electrolyte fuel cell typically operates. 


Pc = —4854.1s7 + 129585, (in Pa) R? = 0.80 
0 < Siw < 0.8 (17) 


It is worthwhile emphasizing that Eq. (17) was utilized solely 
to evaluate the saturation predictions of the traditional Leverett 
approach equipped with J-function for the DM samples tested 
by Ref. [10]. Eq. (17), therefore, should not be considered as 
a generalized relationship appropriate for all type for fuel cell 
DM, although such a relationship is definitely required. 

As expressed, capturing the myriad changes of the transport 
parameters within the pores of the DM is difficult due to the 
minute length scale and the complex heterogeneous nature of the 
DM structure. The distinctive feature of the presented empirical 
fit is that it implicitly accounts for the variations in the rele- 
vant transport properties of the tested DMs including porosity, 
surface tension and contact angle, since it relies on the actual 
experimental capillary pressure data of Ref. [10]. Therefore, 
such a representative curve-fit based on experimental measure- 
ment of the tested DM samples eliminates the requirement of 
porosity, surface tension and J-function as an input. In addi- 
tion the contact angle parameter, which is a required input in 
Leverett approach, is implicitly embedded into the empirical 
curve-fit, thereby accounting for variations in internal contact 
angle caused by the anisotropic nature of the hydrophobic coat- 
ing. This feature enables us to eliminate the need for the selection 
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of a single (and unrealistic) surface contact angle, but Eq. (17) 
should not be considered as a generalized relationship appro- 
priate for all types of fuel cell DM. 


3. Mapping the liquid saturation distribution 
3.1. Modified constitutive relations 


3.1.1. Characteristic relative permeability 

A new characteristic relative permeability correlation appro- 
priate for uncompressed room temperature thin-film fuel cell 
DM was obtained by compiling the experimental capillary 
pressure data presented in Ref. [10] to the empirical models 
described in Section 2.2. The relative permeability correlation 
obtained from these models was fit onto a single curve in the form 
of knw =A(Se-nw)” to represent the overall relative permeabil- 
ity behavior of the tested DMs. Eq. (18) shows the predicted 
relative permeability relation for DMs tested by Ref. [10] at 
room temperature and under no compression. 


kr-nw = (Snw)” É R? = 0.978 (18) 


This relationship should be used with caution, as it is a bulk 
fit for the four different DMs tested by [10], under no compres- 
sion at room temperature. Nevertheless, it is an improvement 
of the existing models applied without linkage to the measured 
experimental data. 
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Fig. 5. (a) Predicted liquid relative permeability curve and (b) percentage devi- 
ation of the Wyllie model from the predicted curve-fit. 


Fig. 5a depicts the predicted relative permeability curve as a 
function of liquid saturation. The predicted relative permeability 
curve for the non-wetting phase shows a modest increase at low 
saturations, because at low saturations only a small portion of 
the pores are occupied by the liquid water, whereas a large por- 
tion of the pores are still available for the gas-phase transport. 
Hence, there is no significant impact of saturation on the reactant 
flow. However, at high saturations, the predicted relative perme- 
ability exhibits a considerable increase with small increases in 
saturation, because as the available pores are nearly all filled by 
the liquid water and gas-phase transport is highly restricted. 

Fig. 5b shows the percentage deviation of the typically used 
Wyllie model (s? model) from the predicted curve-fit at differ- 
ent saturation values. When the present prediction is compared 
with the commonly used Wyllie model (kr-nw = (Spw)°) over 
the entire spectrum of water saturation (0< Snw <1), a signif- 
icant difference between the predicted relative permeability 
values is observed, especially at low water saturation values 
(Snw <0.4). The discrepancy between the predictions dimin- 
ishes at relatively higher saturation values (Spy > 0.6). Physically 
speaking, the coefficients in Wyllie model were generated for 
non-consolidated sands. Thus, the applicability of the Wyllie 
model by itself may not accurately describe the relative perme- 
ability of the fuel cell DM. Currently, the presented approach 
(Eq. (18)) takes into account a broader range of theoretical anal- 
ysis of relative permeability in porous media with the intent of 
providing an improved relative permeability relation and it incor- 
porates direct experimental data obtained in Ref. [10] for various 
common DMs. However, the direct measurement of relative 
permeability for thin fuel cell DM as a function of tempera- 
ture, compression and PTFE content is still needed to verify the 
applicability of present prediction. 


3.2. Implementation of different constitutive relations 


Different constitutive relations were implemented within the 
theoretical model described in Section 2.1 to predict the local 
liquid saturation profiles. The local saturation profiles were pre- 
dicted for two different cases (Case 1 and Case 2 as described 
below), governed by different relations. Furthermore, opera- 
tional conditions and relevant structural properties of thin-film 
porous DM, including net water transport coefficient, current 
density, thickness and hydrophobicity, were examined for these 
two cases in order to capture the significance of the differences 
in liquid saturation profiles. 

In Case 1, the standard Leverett function and Wyllie rela- 
tive permeability relation (kr-nw = s) were used in Eq. (9), 
whereas in Case 2, along with the modified relative permeabil- 
ity (kr-nw = s216), an empirical curve-fit (Eq. (17)) representing 
the measured Pc—s) relation of the DMs tested by Ref. [10] was 
directly integrated into the theoretical model. The relations used 
in each case and the final form of the governing equations for 
these two cases are shown in Table 3. 

Analytical solution of the equations in Table 3 enables deter- 
mination of the local saturation, sı, at any location, x, of the DM. 
The final form of the governing liquid transport equations over- 
all include a number of engineering parameters, including the 
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Table 3 


The constitutive relations used and the final form of the governing equations for Case 1 and Case 2 


Constitutive relations 


Capillary pressure—saturation relationship 


Relative permeability 


1/2 
Case 1 Po=y cos (£) f J(s1) Wyllie model 
Standard Leverett function J(s1) = 1.4175, 2.1205? t 1.2635? knw = sy 
Case 2 Po= —4854.15? + 129585; Present model 
Empirical fit (0< Snw <0.8) R? =0.80 (in Pascals) knw = 52.16 


Final governing equations (Eq. (9)) 


Case 1 


Case 2 


1(2a+1)MWH29 4 
= “\(x) = 0.3545) 


Sog 6 
2Fy cosoV ke 0.8485; j 0.630s7 


H20 
ICat DMW 2 M CX) = 4100.6516 — 2333.754 16 


porosity of the DM, permeability, current density, contact angle, 
surface tension (function of temperature), net water transport 
coefficient and liquid saturation. 


3.3. Validation of the empirical curve-fit 


To evaluate the accuracy of the empirically derived curve- 
fit (Eq. (17)), an additional assessment was performed. 
The saturation prediction of Eq. (17) was compared with 
the saturation profiles predicted by the measured capillary 
pressure-saturation curves of the tested DM samples [10]. The 
capillary pressure—saturation curve of each DM was separately 
implemented into the analytical framework and the correspond- 
ing saturation profiles were calculated at 1 Acm7*, assuming a 
typical DM thickness of 300 wm. Fig. 6 presents the predicted 
saturation profiles of Case 2 by implementing the measured 
capillary pressure—saturation curves of SGL JOBA, SGL 10BB, 
TORAY 090, ETEK-A and the curve-fit given in Eq. (17). The 
saturation profiles over the entire thickness appear to follow sim- 
ilar qualitative trend, predicting very close saturation values. 
More importantly, the curve-fit given in Eq. (17) is observed 
to reasonably well-predict the saturation values for all DM 
samples, with an uncertainly of +6.8%, strongly supporting 
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Fig. 6. Predicted saturation profiles for Case 2 (using empirical curve-fit) imple- 
menting the individual measured capillary pressure curves for all DM samples 
and the curve-fit presented in Eq. (17). 


the effectiveness and applicability of this approach for com- 
paring the experimental data with the prediction of Leverett 
approach. 

In addition, the small deviation between the saturation predic- 
tions indicates that even though the change in PTFE content may 
affect the transport parameters, within the context of the experi- 
mental data provided for the tested DM materials, the predicted 
saturation values appear to be relatively insensitive to the change 
in PTFE within the specified range of DM samples (from 0 to 
5 wt.% PTFE). 


4. Results and discussions 
4.1. Effects of DM properties 


The theoretical model was solved to predict the liquid satura- 
tion profiles for different thicknesses and various hydrophobicity 
levels. The predicted saturation profiles corresponding to dif- 
ferent DM thicknesses were compared. However, the effect of 
hydrophobicity on the predicted local saturations was only per- 
formed for Case 1, since Case 2 is a curve fit for four different 
DMs tailored with similar PTFE content (ranging from 0 to 
5 wt.%) and the contact angle parameter is implicitly embedded 
into this empirical fit, thus not required as an input. The liquid 
saturation was set to zero at the DM-flow channel interface as a 
boundary condition, representing a drainage condition. The con- 
dition of zero liquid surface coverage in the flow channel by liq- 
uid water is reasonably valid for carbon cloth and conditionally 
expected for carbon paper DM under large air stoichometry [6]. 


4.1.1. DM hydrophobicity 

As expected, the hydrophobicity of the DM has a strong 
impact on the liquid saturation distribution predicted by the 
Leverett approach. The liquid saturation profiles for different 
contact angle values were determined for Case | and shown in 
Fig. 7a. At any specified location, the maximum liquid saturation 
is predicted for contact angle 95°, whereas the minimum local 
saturation occurred for contact angle 150°. In addition, for a con- 
stant fractional distance from the catalyst layer, the theoretical 
model was solved for Case 1 to predict the local saturation val- 
ues as a function of contact angle values (shown in Fig. 7b). As 
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Fig. 7. (a) Predicted saturation profiles of Case 1 (using standard Leverett 
approach) for different contact angle values and (b) saturation predictions of 
Case 1 (using standard Leverett approach) at various locations for different 
contact angle values. 


the contact angle increases from 95° to 110°, a relatively severe 
decrease in local saturation values is observed. Above 90°, the 
liquid water transport is enhanced noticeably by an increase in 
hydrophobicity, but with diminishing dependence. Physically, 
the surface adhesion force in carbon fibers is reduced by render- 
ing the DM surface more hydrophobic, causing higher surface 
contact angle. As a result, at higher PTFE loadings of the DM, 
the liquid water molecules on high PTFE loading fibers tend to 
be more unstable, enabling enhanced liquid transport through 
the pores of the DM. 


4.1.2. DM thickness 

Fig. 8 presents the predicted saturation profiles of Case 1 
and 2 for two different DM thicknesses (200 and 300 um). All 
cases (Case | and 2) demonstrate the same qualitative trend, but 
predict lower saturation values for thinner DM at any specified 
location. Physically, for a thinner DM, the liquid water concen- 
tration gradient becomes larger, causing a higher driving force, 
which in turns facilitates the liquid transport. Any decrease in 
DM thickness shortens the path length of the liquid water, thus 
reducing the liquid transport resistance. Therefore, a smaller liq- 
uid saturation difference becomes sufficient to drive the liquid 
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Fig. 8. Predicted saturation profiles for Case 1 (using standard Leverett 
approach) and Case 2 (using empirical curve-fit) for different DM thicknesses. 


water from higher saturation region to lower saturation region 
for a thin DM. 

In terms of saturation predictions of these two Cases (Case 1 
and 2), it can be observed from Fig. 8 that while the cases follow 
the same qualitative trend, the predicted saturation values dif- 
fer significantly as a function of constitutive relation. While the 
capillary-induced transport is defined by Leverett approach in 
Case 1, the discrepancy in saturation predictions between Case 
1 and 2 is noticeable at each location. In Case 1, the standard 
Leverett function (J(s})) is used, however, in Case 2, a polyno- 
mial fit purely derived from experimental data of Ref. [10] has 
been employed. The difference in predicted saturation values 
between each case is relatively more pronounced, even reaching 
a peak value of 0.13 at the MPL interface, as shown in Fig. 8. 

This vast difference in predicted saturation values indicates 
that the standard Leverett function fails to accurately describe 
the capillary transport behavior of the fuel cell DM. This notice- 
able difference can be attributed to the fact that the actual Pc—s) 
relation is directly input into the liquid transport model in Case 2. 
The fuel cell DM has a mixed wettability network, while the Lev- 
erett approach integrates a constant porosity and a representative 
surface contact angle (uniform wettability), which is not the case 
in fuel cell DM. The discrepancy in the model outputs brings 
into question of the appropriateness of assuming constant values 
for porosity and contact angle. An improved approach would 
not only include the rigorous interaction of effective porosity 
and liquid saturation, but would also include the effects of com- 
pressibility, temperature and the mixed wettability (non-uniform 
contact angle) as well. 

Significant structural complexity in the DM leads to diffi- 
culty in defining a meaningful internal contact angle parameter. 
The wettability on any surface can be characterized in terms 
of contact angle, which is closely related to the molecular 
interaction of matters at the contact interface [11,22]. To date, 
empirical studies have been inconclusive regarding the rep- 
resentative average internal contact angle at the microscopic 
level. Due to experimental limitations at the pore-level, exter- 
nal contact angle measurements are commonly employed to 
estimate the internal contact angle, and it is doubtful that an aver- 
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age contact angle approach is appropriate since the connected 
hydrophilic/hydrophobic pores of the DM act in parallel. More- 
over, external contact angle measurements fail to include the 
effects of the different levels of surface energy associated with 
the carbon fibers. Therefore, measured external contact angles 
may not reflect a statistical average of contact angles associated 
within the highly heterogeneous DM pores [23]. 

In addition, the standard Leverett function utilized in Case 
1 accounts for the porous DM in terms of a single dry poros- 
ity value. However, the accumulation of liquid water into the 
available pores and associated change in effective porosity can 
lead to a substantial influence on transport properties of the 
DM [23-25]. Once the available pore space is occupied by the 
liquid water, the liquid flow path resistance increases. There- 
fore, selecting a representative value of DM porosity may fail 
to account for the liquid transport phenomena occurring in fuel 
cells. 


4.2. Operating conditions 


Optimization of the operating conditions is also of great 
importance for the avoidance of flooding problems in the DM. 
The level of liquid saturation in the cathode DM is strongly 
related to the operating current density and net water diffusion 
flux across the membrane [26-29]. The net water transfer coef- 
ficient, a, is strongly coupled with the current density, since the 
rate of water transport across the membrane via diffusion and 
electro-osmotic drag is proportional to the cell operating current 
density. This parameter can be influenced by the selection of DM 
and MPL properties. In the present model, current density and 
net water transfer coefficient are treated as independent param- 
eters to distinguish the significance of each effect separately. 
The theoretical model (Eq. (9)) was solved for different cur- 
rent densities and net water transfer coefficients, while holding 
the other parameters constant. The predicted liquid saturation 
profiles were analyzed and compared. The boundary condition 
is again defined such that the liquid saturation is zero at the 
DM-flow channel interface. 


4.2.1. Net water transfer coefficient 

Fig. 9a shows the effect of net water transfer coefficient, œ, on 
the liquid saturation profile for Case 1. As expected, the present 
model predicts higher liquid saturation in the cathode DM for 
the higher net water transfer coefficient (a = 0.1 versus «œ = 1). As 
the net water transfer rate increases, the cathode side becomes 
more prone to the accumulation of excess liquid water, which in 
turns, increases the local liquid saturation values in the cathode 
DM. 


4.2.2. Operating current density 

Fig. 9b shows the liquid saturation distribution for Case 1 
under different operating current densities (0.2, 1 and 2 Acm7’). 
At any specified location, the predicted saturation value is 
higher for the highest current density (2 Acm~7), whereas the 
lowest liquid saturation is predicted for the lowest current den- 
sity (0.2 A cm7°?), as expected [5]. The same behavior is also 
observed in the MPL. 
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Fig. 9. Predicted liquid saturation profiles for Case 1 (using standard Leverett 
approach); (a) for different net water transfer coefficients, (b) under different 
current densities. 


The saturation profiles predicted by the model for Case 1 
and 2 are compared for a current density of 0.2 A cm7? and a 
value of 0.5 in Fig. 10. As previously discussed, the origin of 
the considerable difference in predicted saturation values arises 
from employing standard Leverett function in Case 1, which 
assumes a constant porosity and contact angle. However, Case 2 
is represented by an experimental Pc—s) curve, which implicitly 
accounts for the variation of DM structural properties such as 
internal contact angle and effective porosity. This remarkable 
deviation in the predicted local saturation values between Case 
1 and 2 further demonstrates the inapplicability of using the 
standard Leverett approach. 


4.3. Case sensitivity analysis 


A sensitivity analysis was performed to ascertain how each 
specified case prediction depends upon the DM thickness, the net 
water transfer coefficient and the current density, and to assess 
the range over which the model predictions become sensitive to 
changes in these parameters. 

As a first step, the functional dependence of the liquid satura- 
tion as a function of a chosen input parameter was determined at 
a specified location (e.g. a fractional distance from the catalyst 
layer of 0.8). Then the derivative of this function with respect to 
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Fig. 10. Predicted saturation profiles for Case 1 (using standard Leverett 


approach) and Case 2 (using empirical curve-fit) for 200 uym DM and a@ value 
of 0.5; (a) at current density 1.0A cm~?, (b) at current density 0.2 A cm. 


the chosen input parameter was computed to evaluate the sen- 
sitivity of the model prediction corresponding to any change in 
that input parameter. A greater slope implies a greater sensitiv- 
ity, as shown in Eq. (19). The derivatives with respect to the 
specified inputs are also compared between Case 1 and Case 2. 


(S1)Case-j = Func(X) 


ORO 
ax Case- j ax Case-i 


=> Case-j is more sensitive than Case-i. (19) 


at a specified location 


The first investigated input parameter is the DM thickness. 
Fig. 11a depicts the predicted saturation profiles for two dif- 
ferent DM thicknesses. As the DM thickness increases from 
200 to 300 um, the predicted saturation values at each location 
increases for all cases. It is observed that the average percent- 
age increase in predicted saturation values as the thickness is 
increased from 200 to 300 um is found to be 10.3% for Case 1 
and 10.5% for Case 2, approximately equal. The sensitivity of 
the case predictions corresponding to any change in DM thick- 
ness was also investigated at a fractional distance from catalyst 
layer of 0.8. Fig. 11b shows the rate of change of the predicted 
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Fig. 11. (a) Predicted saturation profiles for Case 1 (using standard Leverett 
approach) and Case 2 (using empirical curve-fit) for different DM thicknesses, 
(b) rate of change of saturation predictions for the two cases with varying DM 
thickness. 


saturation values corresponding to the DM thicknesses from 100 
to 500 um. As seen from the graph, the magnitudes of the slopes 
are on the order of 1074, meaning that DM thickness seems to 
have minor effect on case predictions. 

A case sensitivity analysis of the net water transfer coefficient 
was also performed. Two different values of net water transfer 
coefficient (œ = 0.1 and 1.0) were implemented and the predicted 
saturation profiles are shown for Case 1 and 2 in Fig. 12a. The 
increase in net water transfer coefficient from 0.1 to 1 appears to 
cause a slight increase in the predicted saturation values. How- 
ever, the amount of increase is different for Case 1 and 2. In order 
to identify the relative significance of changing the net water 
transfer coefficient on case predictions, the sensitivity analysis 
was performed at a specified location (a fractional distance from 
the catalyst layer of 0.8). Fig. 12b depicts that the sensitivity of 
the each case is reduced as the net water transfer coefficient 
increases from 0.1 to 0.5 and then starts to increase as the net 
water transfer coefficient increases from 0.5 to 1.0. As seen from 
Fig. 12b, Case 2 exhibits relatively smaller slope change, rang- 
ing from 0.013 to 0.085, whereas Case 1 has higher slope values 
(varies from 0.03 to 0.2) over the entire net water transfer coef- 
ficient, meaning that model predictions for Case 1 is relatively 
more sensitive, (by an average factor of 2.2) than that of Case 2. 
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Fig. 12. (a) Predicted saturation profiles using Case 1 (using standard Leverett 
approach) and Case 2 (using empirical curve-fit) for different net water transfer 
coefficients, (b) rate of change of saturation predictions for all cases with varying 
net water transfer coefficients. 


The final investigated input parameter is current density. Two 
different operating current densities ([=0.2 and 1.0 A cm~?) 
were implemented in two cases and the behavior of the model 
predictions were analyzed. Fig. 13a shows the predicted liquid 
saturation distribution for Case 1 and 2 under different current 
densities. As expected, the increase in current density gives 
rise to the predicted saturation values at each location due to 
the enhanced water generation. However, when the amount of 
average increase in predicted saturations is compared, Case 1 
exhibits relatively smaller percentage increase (59%) compared 
to Case 2 (65%), in response to the change in current density. 
Fig. 13b depicts the rate of change of the predicted saturation 
values as a function of current density at a specified location 
(a fractional distance from catalyst layer of 0.8). As seen from 
Fig. 13b, Case 1 is observed to be more sensitive compared to 
Case 2 with an average factor of 2.08, corresponding to any 
changes in current density over the range of J=0-2 Acm7?. 
The sensitivity is more pronounced for these two cases when 
the current density varies from 0 to 0.5 A cm~? and from 1.5 
to 2Acm~?. Conversely, the predicted sensitivity for all cases 
becomes relatively small when the current density is in the range 
of 0.5<I< 1.5 A cm™°?. 

The magnitude of the slope with respect to any specified 
input parameter can be used to compare the relative significance 
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Fig. 13. (a) Predicted saturation profiles using Case 1 (using standard Leverett 
approach) and Case 2 (using empirical curve-fit) under different current density, 
(b) rate of change of saturation predictions for all cases under various operating 
current densities. 


of each input parameters on the model prediction. Among the 
input parameters of concern (thickness, current density and net 
water transfer coefficient), the current density seems to be the 
most crucial input parameter for both Case 1 and Case 2, when 
the magnitude of the rate of change is compared between the each 
input parameters. Based on these analyses, it can be concluded 
that the model prediction for all cases is relatively more sensitive 
to any change in current density rather than the changes in the 
DM thickness and net water transfer coefficient. 


5. Conclusions 


The effectiveness of the traditional Leverett approach 
(equipped with J-function) to describe the capillary-induced 
flow in thin-film fuel cell DM was investigated. An 
empirical polynomial fit (Pc-sı) describing the capillary 
pressure-saturation relation of the DMs listed in Table 2 was 
derived based on the experimental capillary pressure—saturation 
curves provided by Ref. [10]. In addition, a new representa- 
tive relative permeability correlation (kr—nw = s216) compiling 
a broader range of theoretical basis for the tested fuel cell media 
by Ref. [10] was obtained. These correlations, along with the 
Leverett approach equipped with J-function, were integrated 
into an analytical framework to predict liquid saturation pro- 
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files for differently engineered DMs under various operating 
conditions. The variations in liquid saturation predictions asso- 
ciated with employing different correlations were assessed and 
the effectiveness of the traditionally used Leverett function was 
determined to be poor. 

The Leverett approach consistently over-predicts the DM 
saturation values compared to the ones obtained with the empir- 
ical fit of the experimental data. The significant variation in 
the predicted saturation profiles suggests that using the tradi- 
tional Leverett approach is indeed inadequate to describe the 
capillary transport characteristic of a fuel cell DM, indicating 
the need to develop a modified approach appropriate for thin- 
film fuel cell media with mixed wettability. Finally, a sensitivity 
analysis was performed to ascertain how each specified corre- 
lations prediction depends upon the different input parameters. 
The model predictions are found to be relatively more sensitive 
to any change in current density rather than the changes in the 
DM thickness and net water transfer coefficient over the range 
of values tested. 

Note that the presented correlations were generated based 
on the data [10] for a set of DMs tailored with PTFE ranging 
from 0 to 5 wt.% at room temperature under no compression. 
However, hydrophobicity of the DM, along with the operational 
environments such as temperature and compression needs to be 
investigated in a broader range to further isolate the effects of 
these parameters on the capillary transport characteristics of a 
fuel cell DM. 
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